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ABSTRACT 

We use high-resolution three-dimensional adaptive mesh refinement simulations to investigate the 
interaction of high-redshift galaxy outfiows with low-mass virialized clouds of primordial compo- 
sition. While atomic cooling allows star formation in objects with virial temperatures above 10^ 
K, "minihaloes" below this threshold are generally unable to form stars by themselves. However, 
these objects are highly susceptible to triggered star formation, induced by outfiows from neighboring 
high-redshift starburst galaxies. Here we conduct a study of these interactions, focusing on cooling 
through non-equilibrium molecular hydrogen (H2) and hydrogen deuteride (HD) formation. Tracking 
the non-equilibrium chemistry and cooling of 14 species and including the presence of a dissociating 
background, we show that shock interactions can transform minihaloes into extremely compact clus- 
ters of coeval stars. Furthermore, these clusters are all less than ^ lO^M©, and they are ejected from 
their parent dark matter halos: properties that are remarkably similar to those of the old population 
of globular clusters. 

Subject headings: galaxies: formation- galaxies: star clusters - globular clusters: general - astrochem- 
istry - galaxies: high-redshift - shock waves 



1. INTRODUCTION 

A generic prediction of the cold dark matter model is 
a large high-redshift population of gravitationally-bound 
clouds that are unable to form stars. Because atomic 
H and He line cooling is only effective at temperatures 
above 10^ K, clouds of gas and dark matter with virial 
temperatures below this threshold must radiate energy 
through dust and molecular line emission. While the 
levels of H2 left over from recombination are sufficient 
to cool gas in the earliest structures {e.g. Abel et al. 
2002; Bromm et al. 2002), the resulting 11.20 - 13.6 eV 
background emission from the stars in these objects {e.g. 
Haiman et al. 1997; Ciardi et al. 2000; Sokasian et al. 
2004; O'Shea & Norman 2007) is likely to have quickly 
dissociated these trace levels of primordial molecules 
(Galli & Palla 1998). And although an early X-ray back- 
ground could have provided enough free electrons to pro- 
mote H2 formation, the relative strength between these 
two backgrounds is uncertain, and it is unlikely that 
the background was strong enough to balance ultravi- 
olet (UV) photodissociation. Even if there were some 
trace amount of H2 in these clouds, it is likely to be in 
such a small abundance as to not impact their structure 
(Whalen et al. 2008a; Ahn et al. 2009). 

The result is a large number of dark matter halos that 
were massive enough to overcome the thermal pressure of 
the primordial intergalactic medium and retain their gas, 
but not massive enough to excite the radiative transitions 
necessary to cool this gas into stars. These "minihalos," 
whose virial temperatures were Tvir < 10^ K and whose 

total masses were between 10^ and 10^"^ M© would have 
remained sterile until some outside force either disrupted 
them, or more interestingly, disturbed them so as to cat- 
alyze coolant formation. In fact, two possible coolant 
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formation methods have been considered in detail in the 
literature: ionization fronts and shock fronts. 

In the case of ionization fronts, such as would occur 
during the epoch of reionization, high-energy photons 
emitted from galaxies or quasars interact with the neu- 
tral atomic minihalo gas. Bond, Szalay & Silk (1988) 
originally discussed how the resulting photoionization 
would expel the gas contained in a minihalo by sud- 
denly heating it to T ^ 10^ K, as would be the case 
in the optically-thin limit. On the other hand, Cen 
(2001), used simple analytic estimates to argue that ion- 
ization fronts would cause non-equilibrium H2 formation 
and the collapse of the gas inside the gravitational po- 
tential. However, Barkana & Loeb (1999) studied mini- 
halo evaporation using static models of uniformly illu- 
minated spherical clouds, accounting for optical depth 
and self-shielding effects, and showed that the cosmic UV 
background boiled most of the gas out of these objects. 
Later, Haiman, Abel, & Madau (2001) carried out three- 
dimensional (3D) hydrodynamic simulations assuming 
the minihalo gas was spontaneously heated to 10^ K, 
also finding quick disruption. Finally, full radiation- 
hydrodynamical simulations of ionization front-minihalo 
interactions were carried out in Iliev et al. (2005) and 
Shapiro et al. (2004; see also Shapiro, Raga & Mellema 
1997; 1998). These demonstrated that intergalactic ion- 
ization fronts decelerated when they encountered the 
dense, neutral gas inside minihaloes and were thereby 
transformed into D-type fronts, preceded by shocks that 
completely phot oe vapor ated the minihalo gas. 

A second and more promising avenue for coolant for- 
mation is the interaction between galactic outfiows and 
minihaloes. These galaxy-scale winds, which are driven 
by core-collapse supernova and winds from massive stars, 
are commonly observed around dwarf and massive star- 
bursting galaxies at both low and high redshifts {e.g. 
Lehnert & Heckman 1996; Franx et al. 1997; Pettini et al. 



1998; Martin 1999; 1998; Heckman et al 2000; Veilleux et 
al. 2005; Rupke et al. 2005), and a variety of theoretical 
arguments suggest that these galaxies represent only the 
tail end of a larger population of smaller "pre-galactic," 
starbursts that formed before reionization (Scannapieco, 
Ferrara & Madau 2002; Thacker, Scannapieco, & Davis 
2002). Furthermore, the interstellar gas swept up in a 
starburst-driven wind can effectively trap the ionizing 
photons behind it (Fujita et al. 2003), meaning that at 
high redshifts, many intergalactic regions may have been 
impacted by outflows well before they were ionized. 

Such ^ 100—300 km/s shocks can cause intense cooling 
through two mechanisms: (i) the mixing of metals with 
ionization potentials below 13.6 eV (Dalgarno & McCray 
1972), which allow for atomic line cooling even at tem- 
peratures below 10^ K; and (ii) the formation of H2 and 
HD by nonequilibrium processes (Mac Low & Shull 1986; 
Shapiro & Kang 1987; Kang et al. 1990; Ferrara 1998; 
Uehara & Inutsuka 2000), which allow for molecular line 
cooling associated vibrational and rotational transitions 
(Palla & Zinnecker 1988). In fact, Scannapieco et al. 
(2004) showed that these effects were so large that shock 
interactions could induce intense bursts of cooling and 
collapse in previously "sterile" minihalo gas. Using sim- 
ple analytic models, they found that the most likely re- 
sult was the formation of compact clusters of coeval stars, 
although they also emphasized the importance of multi- 
dimensional numerical studies to confirm this result. 

In this paper, we present the first in a series of 3D nu- 
merical studies to capture the physical interactions be- 
tween primordial minihalos and supernova- driven galac- 
tic outflows. While triggered star formation has been 
simulated in low-redshift intergalactic clouds impacted 
by radio jets {e.g. van Breugel et al. 1985; Fragile et al. 
2004; Wiita 2004; Klamer et al. 20004), this has never 
been numerically simulated in the high-redshift, mini- 
halo case. In this paper, we focus on the physics of non- 
equilibrium II2 and HD chemistry and associated cooling 
in determining the properties of the post-shock minihalo 
gas. In particular, we show that molecule formation and 
the resulting cooling is strong enough to induce rapid 
minihalo collapse and star formation, leading to compact 
stellar clusters. 

The structure of this paper is as follows. In §2 we out- 
line the physical components of the model, focusing on 
our primordial H2 and HD chemical network and cool- 
ing routines and their respective tests. In §3 we outline 
the general model used for the galactic outflow and the 
minihalo, and in §4 we present the simulation results and 
relate the resulting stellar clusters to local observations. 
Conclusions are given in §5. 

2. NUMERICAL METHOD 

All simulations were performed with FLASH version 
3.1, a multidimensional adaptive mesh refinement hydro- 
dynamics code (Fryxell et al. 2000) that solves the Rie- 
mann problem on a Cartesian grid using a directionally- 
split Piecewise-Parabolic Method (PPM) solver (Colella 
& Woodward 1984; Coleha & Glaz 1985; Fryxeh, Miiller, 
& Arnett 1989). Furthermore, unlike earlier versions of 
the code, FLASH3 includes an effective parallel multi- 
grid gravity solver as described in Ricker (2008). How- 
ever, before shock-minihalo interactions could be simu- 
lated, two capabilities needed to be added: nonequili- 



birum primordial chemistry, and cooling from atoms and 
from molecules produced in these interactions. In this 
section we describe our numerical implementation of each 
of these processes, along with the tests we carried out be- 
fore applying the code to shock-minihalo interactions. 

2.1. Chemistry 

As the minihalos we are considering in this paper are 
made up of primordial gas, their chemical makeup is 
highly restricted, with contributions from only hydro- 
gen, helium, and low levels of deuterium. Yet even these 
three isotopes can exist in a variety of ionization states 
and molecules and are thus associated with a substan- 
tial network of chemical reactions that must be tracked 
throughout our simulations. 

2.1.1. Implementation 

The chemical network that was implemented into 
FLASH is outlined by Glover & Abel (2008, hereafter 
GA08). Throughout our simulations we track three 
states of atomic hydrogen (H, H+, & H~) and atomic 
deuterium (D, D+, & D~), three states of atomic helium 
(He, He+, &, He++), two states of molecular hydrogen 
(H2 & H2") and molecular hydrogen deuteride (HD & 
HD+), and electrons (e~). For simplification, any re- 
action that involved molecular deuterium (D2) and all 
three-body reactions were neglected. As stated in GA08, 
the very small amount of D2 and D2" produced makes any 
cooling by these molecules irrelevant, while three-body 
reactions only become important at n ^ 10^ cm~^ (e.^. 
Palla et al. 1983), many orders of magnitude denser than 
the conditions considered here. With these constraints, 
a total of 84 reactions were used out of the 115 described 
in GA08. 

Photodissociation rates due to an external radiation 
field were also included as given in Glover & Savin (2009). 
These rates are calculated assuming a Teff = 10^ K black- 
body source and their strength is quantified by the flux 
at the Lyman limit, J(z^a) = 10~^^J2i erg s~^ cm~^ 
Hz~^ sr~^. Note that once H2 and HD are produced 
in sufficient quantities, some molecules are self-shielded 
from the background radiation. However, for simplic- 
ity, we considered only the case where there was no self- 
shielding, and thus our results place an upper limit on the 
effect of a dissociating background. This process adds an 
additional 7 reactions for a total of 91 reactions in the 
chemical network. 

In reactions that involve free electrons recombining 
with ions, there are two possible choices for the reaction 
rate, depending on the overall optical depth of the cloud 
to ionizing radiation. In the optically-thin case (Case 
A; Osterbrock 1989) ionizing photons emitted during re- 
combination are lost to the system, while in the optically- 
thick case (Case B), ionizing photons are reabsorbed by 
neighboring neutral atoms, which have the effect of low- 
ering the recombination rates by essentially not allowing 
recombination to the ground state. There are three re- 
actions (H+ + e~ ^ H + 7, He^ + e~ ^ He + 7, and 
He^^ + e~ ^ He^ + 7) where this is a concern, and as 
we shall see below, in all cases our clouds are optically 
thick, such that Case B recombination rates are appro- 
priate. 
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The binding energy from each species is also important 
in the total energy budget and on the evolution of the 
gas. In FLASH these are defined such that all of the 
neutral species (H, D, e~, & H) have binding energies 
equal to zero. As the gas is heated and the atomic species 
begin to ionize, the endothermic reactions remove the 
binding energy between the nucleus and electron (s) from 
the internal energy of the gas. For H and D this requires 
13.6 eV, while the ionized states for Helium (He+ and 
He++) have ionization potentals of 24.5 eV and 79.0 eV 
respectively. H~ and D~ are only weakly bound and 
have similar binding energies of 0.75 eV. Finally H2 and 
HD have binding energies of 4.4 eV, and H2" and HD+ 
have binding energies of 10.9 eV, somewhat lower than 
the atomic species. 

To describe the evolution of our 14 species, we enu- 
merate them with an index i such that each has Zi pro- 
tons and Ai nucleons, following the structure and syntax 
from Timmes (1999). Next we consider a gas with a to- 
tal mass density p and temperature T and denote the 
number and mass densities of the ith isotope as n^ and 
Pi, respectively. For each species we also define a mass 
fraction 

Xi = pi/p = niAi/{pNA), (1) 

where Na is Avogadro's number, and we define the molar 
abundance of the ith species as 



Y,=XjA, = nJ{pNA), 



(2) 



^N 



where conservation of mass is given by ^- Xi = 1. Each 
of the 14 species can then be cast as a continuity equation 
in the form 

^^f = R. (3) 

where Ri is the total reaction rate due to all the binary 
reactions of the form i -\- j -^ k -\- 1^ defined as 



Ri 



Y,YiYk\kj{l) 



YiYjXjk{i), 



(4) 



where Xkj and Xjk are the creation and destruction chem- 
ical reaction rates for a given species. If the species in 
question is affected by UV background radiation, the 
continuity equation takes the following form. 



R, = y2YiYkXkj{l) 



j,k 



YiYjXjkii) - YiJii^a), (5) 



where the last term accounts for the amount of these 
species that are destroyed by the background radiation, 
J(z^c^). Throughout our simulations, changes in the num- 
ber of free elections are not calculated directly, but rather 
at the end of each cycle we use charge conservation to 
calculate their molar fraction, as 

Yelec = ni++>D++>^HD++^H++^He++2yHe++->H-->D- 

(6) 

Because of the often complex ways that the chemi- 
cal reaction rates depend on temperature and the intrin- 
sic order of magnitude spread in the rates, the resulting 
equations are 'stiff,' meaning that the ratio of the min- 
imum and maximum eigenvalue of the Jacobian matrix, 
Jij = dYi/dYj^ is large and imaginary. This means that 



implicit or semi-implicit methods are necessary to effi- 
ciently follow their evolution. To address this problem, 
we arrange the molar fractions of the 13 species, exclud- 
ing e~, into a vector Y, and solve the resulting system 
of equations using a 4^^ order accurate Kaps-Rentrop, 
or Rosenbrock method (Kaps & Rentrop 1979). In this 
method, the network is advanced over a time step h via 



Y"+i=Y" + ^6iAi, 



(7) 



where the Ai vectors are found by successively solving 
the four matrix equations 

(iM-J)-Ai=/(Y-), (8) 

(iM -J) '^2 = /(Y^ + ^21 Ai) + C21 A1//1, (9) 
{i/jh - J) • A3 = /(Y- + aai Ai + a32 A2) 

+ (C3i Ai + C32 A2)//i, and (10) 

(i/7/1 - J) • A4 = /(Y^ + a4i Ai + a42 A2 + a43 A3) 

+ (C41 Ai + C42 A2 + C43^3)/h. (11) 

Here bi^j^aij, and Cij are fixed constants of the 

method, /(Y) = Y, 1 is the identity matrix, and J is 
the Jacobian matrix. Note that the four matrix equa- 
tions represent a staged set of linear equations and that 
the four right hand sides are not known in advance. At 
each step, an error estimate is given for the difference be- 
tween the third and fourth order solutions. For compari- 
son we also carried out tests, using a multi-order Bader- 
Deuflhard method (Bader & Deuflhard 1983). However 
in the end, the Rosenbrock method was chosen over this 
method because of its efficiency and speed. 

As the species evolve, the temperature of the gas 
changes from the release of internal energy from recombi- 
nations or the loss of internal energy from ionizations and 
dissociations. These changes can in turn affect the reac- 
tion rates. Thus to ensure the stability of the chemistry 
routine while at the same time allowing the simulation 
to proceed at the hydrodynamic time-step, we developed 
a method of cycling over multiple Kaps-Rentrop time 
steps within a single hydrodynamic time step. Here we 
estimated an initial chemical time step of each species as 



'^cherQji ^chem 



Yi+ 0.1YH+ 
Y 



(12) 



where achem is a constant determined at runtime that 
controls the desired fractional change of the fastest evolv- 
ing species. The change in molar abundances, F^'s, were 
calculated from the ordinary differential equations that 
make up the chemical network, and the molar fractions 
of each species Y^'s are given by the current values. 
In both the tests and simulations, we chose a value of 

Q^chem = 0.5. 

Note that we offset the subcycling time step by adding 
a small fraction of the ionized hydrogen abundance to eq. 
( 12 ). This is because there are conditions where a species 



is very low in abundance but changing very quickly, for 
example, rapid ionization of atomic species, which will 
cause the subcycling to run away with extremely small 
time steps. In regions in which most species are neutral, 
this has little effect since the chemical time step is likely 



longer than the hydrodynamical one, and in regions in 
which abundances are rapidly changing, then this ex- 
tra term buffers against very small times steps. It also 
prevents rapid changes in internal energy as energy is re- 
moved as atomic species are ionized and gained as they 
recombine. 

Once calculated, these species time steps are compared 
to each other and the smallest time step, associated with 
the fastest evolving species, is chosen as the subcycle time 
step. If this is longer than the hydrodynamic time step, 
the hydrodynamic time step is used instead and no ad- 
ditional subcycling is done. If subcycling is required, the 
species time step is subtracted from the total hydrody- 
namic time step and the network is then updated over the 
chemical time step. The species time steps are recalcu- 
lated after each subcycle and compared to the remaining 
hydrodynamic time step. This is repeated until a full hy- 
drodynamic time step is completed, as is schematically 
shown in Fig. [l] 
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and the network can be bypassed. If the temperature is 
between 2.0 x 10^ and 1.0 x 10^ K , then we 'prime' the 
solutions and ionize 5% of available neutral hydrogen, 
5% of neutral helium (4.5% into singly ionized helium 
and 0.5% into doubly ionized helium), before entering 
the iterative solver, to help accelerate the routine to- 
wards the correct solution. Finally if the temperature is 
less than 50 K, then all species are kept the same, and 
no reactions are calculated. This is done because cool- 
ing and chemistry rates become unimportant at such low 
temperatures. It also has the benefit of speeding up the 
simulation slightly as very little time is spent in either 
the cooling or chemistry routines. In all other cases, the 
full network is evolved without alterations. 

2.1.2. Chemistry Tests 




Time (s) 



Time (s) 



Time (s) 



Fig. 2. — Chemical evolution tests. Column 1 shows the T = 
10^ K case, column 2 shows the T = 10^ K case, and column 
3 shows the T = 10^ K. Time is given on the x-axis and the 
number density of each species divided by the total number density 
of hydrogen is given on the y-axis. The blue lines correspond to 
the n = 0.01 cm~^ case, red to the n = 0.1 cm~^ case, green to 
the n = 1.0 cm~^ case, magenta to the n = 10.0 cm~^ case, and 
teal to the n = 100.0 cm~^ case. The solid lines are results from 
FLASH and the dashed lines are results from G09. 



Fig. 1. — Schematic view of chemistry subcvcli ng. First, the 
chemical time step, r{c), is calculated using eqn. ( |12| >. If this is 
larger than the hydrodynamic time step, then the evolution time 
step, r(e), is set to the hydrodynamic time step. Else, r(e) is set to 
the chemical time step. The network is then evolved for r(e) and 
the remaining time step t(/z) is calculated. If this is zero then we 
proceed to the next step, else we cycle back through the network 
with the updated abundances. This loop continues until the full 
hydrodynamic time step is covered. Note that after every chemical 
network iteration, the cooling routine is called. 



In cases in which the gas is extremely hot or cold, the 
chemical make-up can be determined directly from the 
temperature, avoiding the need for matrix inversions. If 
the temperature is above 10^ K then all atomic species 
become ionized and all molecular species are dissociated. 



To test our implemented chemical network, we carried 
out a series of runs in which initially dissociated and 
ionized gas was held at constant temperature and den- 
sity for 10^^ seconds. A small initial time step {to ~ 
10^ s) was used and allowed to increase up a maximum 
time step of 10^^ seconds. Models were run with to- 
tal hydrogen number densities varying from 0.01 cm~^ 
to 100 cm~^, and temperatures ranging between 10^ K 
and 10^ K, and no external radiation. In each case, the 
results were compared to the results of a different im- 
plementation of the same chemical network within the 
ENZO code (Glover 2009, private communication, G09), 
yielding the molar fractions shown in Figure |2j In this 
figure, the three columns correspond to runs with differ- 
ent temperatures, the curves corresponds to runs with 
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different densities, and the rows correspond to the evo- 
lution of different species. 

The match between our tests and the numerical results 
from G09 is excellent. In all cases and at all tempera- 
tures, the curves closely track each other, in most cases 
leading to curves that are indistinguishable. Although 
the abundances of several species change by many or- 
ders of magnitude throughout the runs, the two methods 
track each other within to 10% in all cases except for H2 
at lO^K, which is unimportant as a coolant but never- 
theless consistent within a factor of 1.5 at all times. Fur- 
thermore, this agreement between methods is also seen 
for deuterium species, which are not shown in this fig- 
ure as they follow H exactly, maintaining a g^ ratio 
between both species at all times. 

At T = 100 K, all ionized species quickly recombine 
with the free electrons to form neutral atoms. However, 
even during this relatively quick transition from ionized 
to neutral, H+ and H~ ions (not shown) persist for long 
enough to catalyze the formation of substantial amounts 
of molecular gas, leading to final H2 molar fractions of 
^ 10~^. At T = 1000 K, the evolution is very similar to 
the T = 100 K case, although the species do not reach 
equilibrium as quickly, leading to even higher levels of H2 
formation. Finally, at T = 10^ K, it takes even longer for 
the ionized hydrogen to recombine, but in this case, less 
molecular species are formed, as collisional dissociation 
of H2 and HD are more prevalent, limiting the maximum 
amount of these species. 

Also apparent in these plot is the dependence of the 
species evolution on the density of the gas. Chemical 
reactions are fundamentally collisional processes whose 
rates are quadratic in number density. Thus, as we are 
not considering three-body interactions, the timescale as- 
sociated with chemistry should decrease linearly with the 
density. This is seen for all temperatures and species 
shown in Figure [2| as in every case each line is sepa- 
rated from its neighbor by a factor of 10 in time, exactly 
corresponding to the density shift between cases. 

2.1.3. Effect of the Background Radiation 

Background radiation with photon energies between 
11.2 and 13.6 eV can excite and dissociate molecular hy- 
drogen. In the absence of other coolants, this can have 
drastic effect on the evolution of the cloud. Two extremes 
are immediately apparent, a strong background case in 
which any II2 or HD formed is quickly dissociated, and a 
background-free case in which no molecules are photodis- 
sociated. A simple test was constructed to study the ef- 
fect of the background and determine a fiducial value for 
J21. J21 was varied between and 1 at five different val- 
ues. For each value of J21 the number density was varied 
between n = 10~^ and 1.0 cm~^. Each test was run at 
a constant temperature and constant density with evolv- 
ing chemistry and no cooling. The results are given in 
Figure [3] From this, we determine that only background 
levels at or above J21 = 0.1 give an appreciable differ- 
ence in the abundance of H2 and HD over a megayear 
timescale, which as we shall see below, is the timescale 
of shock-minihalo interactions. At the same time, J21 
= 0.1 provides a reasonable upper limit to the level of 
background expected before reionization {e.g. Ciardi & 
Ferrara 2005). Therefore, we use this value as a fiducial 
value in the simulations with a background. 



n = 0.1 cm-3 




1 1.5 2 

Log[Myr] 



Fig. 3. — Comparison of different UV backgrounds. The dotted 
line is the comparison from G09, the sohd hne is J21 = 0, the 
short-and-long-dashed hne is J21 = 10~^, the dot-long-dahsed hne 



is J21 



10- 



the dot-short-dashed hne is J21 = 10 , the long- 



dashed line is J21 = 10"-*^, and the short-dashed line is J21 = 1.0. 
Time is given on the x-axis and number density of each species 
normalized by the number density of neutral hydrogen is given on 
the y-axis. Note that the solid line and dotted lines coincide with 
each other, demonstrating that we recover the expected results in 
the background- free case. 

2.2. Cooling 

The second major process added to the code was ra- 
diative cooling, which was divided into two temperature 
regimes. At temperatures > lO^K, cooling results mostly 
from atomic lines of H and He, with bremmstrahlung ra- 
diation also becoming important at temperatures above 
lO^i^. Below 10^ K, on the other hand, the net cooling 
rate is determined by molecular line cooling from H2 and 
HD, which, as it is an asymmetric molecule, can radiate 
much more efficiently than H2, and thus can be almost 
as important although it is much less abundant. Cooling 
from H2 operates down to T < 200 K and to number den- 
sities n > 10"^ cm-3 (Glover & Abel 2008; Gahi & Paha; 
1998), while HD which can cool the gas to slightly lower 
temperatures and to higher number densities (Bromm, 
Coppi, & Larson 2002). As we are restricting ourselves 
to primordial gas in this study at any given temperature 
the overall cooling rate, Axotah is the combination from 
both regimes. 



Axotal — AAtomic + AMolecular- 

Each cooling rate has the form: 



A., 



/ vo I v n y^o n ^ 



(13) 



(14) 



where Aij is the energy loss per volume due to species i 
and j. Hi and Uj are the number densities of each species, 
and Xij is the cooling rate in ergs cm~^ s~^. Cooling 
rates for the collisional excitation between H2 and H, H2, 
H+, and e~ and between H2" and H or e~ are taken from 



s 




Ei = Ei-s-T(e) 
Calc:T 



yf 


Ei = (l-acoo,)Ei 

T = (l-aco=,)T 

T(e) = T(e)-T(cl) 



Done. 
Back to Chemistry. 



Fig. 4. — Schematic view of the coohng subcycle. The time over 
which chemistry evolves r(e) is used as the initial time step. T his i s 
compared against t{cI) the cooling time step, as given by eq. ( |15| . 
If the cooling time step is shorter than the evolved time step, then a 
portion of the internal energy and temperature are subtracted and 
the evolved time step is updated. The cooling time step is then 
recalculated. If the cooling time step is longer than the evolving 
time step then the internal energy is directly updated and used 
to calculate the new temperature. Once this is done, cooling is 
complete and we return to the chemistry routine 

GA08. The cooling rate for the cohisional excitation be- 
tween HD and H is taken from Lipovka, Nunez-Lopez, & 
Avila- Reese (2005). Finahy, coohng rates from Hydrogen 
and Hehum atomic hnes are calculated using CLOUDY 
(Ferland, G.J., et al 1998). In calculating these rates, we 
followed the procedure described in Smith et al. (2008) 
and used the "coronal equilibrium" command which con- 
siders only cohisional ionization. The cooling curve was 
calculated assuming case B recombination for the recom- 
bination lines of hydrogen and helium, as discussed fur- 
ther in §3.1. 

Any cooling routine contains a natural timescale that 
relates the total internal energy to the energy loss per 
time: 

acooi X Ei 

^cool = -. , (15) 

s 

where acooi is a constant between and 1, in all cases set 
at 0.1, Ei is the internal energy, and s is the energy loss 
per time. Cooling rates are very dependent on temper- 
ature and species abundances and these quantities can 
change rapidly over a single chemical time step. 

A method of subcycling over cooling time steps was 
developed to ensure that the correct cooling rates are 
used. An initial coolin g ti me scale is calculated assuming 
<^cooi = 0-1 using eqn. (15) which is then compared to the 
chemical time step. If TcooI is smaller than the fraction 
of the chemistry time step then that fraction of energy 
is subtracted from the internal energy and temperature. 
The cooling rate and cooling time step is recalculated 
with the updated temperature. This continues until the 
chemistry time step is reached. This is schematically 
given in Fig. [4j 
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Fig. 5. — Coohng tests. The sohd hnes are taken from Prieto et 
al. (2008) and compared to our model. The blue curves correspond 
to a number density n = 1.0 cm~^, red to n = 10.0 cm~^, and 
green to n = 100.0 cm"*^. The temperature is not allowed to go 
below 50 K. 



As a test of our cooling routines, we reproduced the 
example curves given in Prieto et al. (2008). In this 
work, the authors present the effects of H2 and HD cool- 
ing in a primordial gas. The gas begins at an initial 
temperature of T = 500 K with initial number densi- 



ties, relative to hydrogen: nH+ = 10 ; nH 



-5., 



n„+ 



10 ^;nHD 



2.2.1. Cooling Tests 



10" ;nH2 = 10"'^;nD = 10"";nD+ 
10~^;nHD+ = 10~^^;^He+ = ^He++ = 0-Oi ^^^ with ini- 
tial hydrogen and helium densities of pn = 0.75 x ptot 
and pHe = 0.24 x ptot^ where ptot is the total baryonic 
matter density. 

Three models were run with total number densities of 
Utot = 1.0,10.0, and 100.0 cm~^. Cooling was tracked 
for 10^ yrs with chemistry evolving simultaneously. The 
results of this calculation are shown in Figure |5J which 
indicates good agreement with Prieto et al. (2008). It 
should be noted that temperature evolution in this plot 
has a linear dependence of the number density of the gas. 
For example, a gas with ten times the number densities 
of another gas will cool ten times quicker. This is again 
because most of the cooling is coming from the collisions 
between two species, in this case H2 or HD and H. 

As mentioned above, HD can be more important than 
H2 for gas cooling at higher densities and colder tem- 
peratures. To determine whether or not HD cooling is 
important in this simulation, we apply the cooling test 
to two different scenarios. First, we use the same ini- 
tial abundances as described above and second, using 
primordial abundances with a small fraction (0.01%) of 
each atomic species ionized. Each test was run twice, 
once with deuterium and once without. The results of 
these tests is given in Figure [6J At high number densi- 
ties, HD cooling does not have a perceivable effect. At 
intermediate temperatures, HD cooling is important for 
a gas with the initial abundances from the cooling tests 
Finally, at low temperatures, HD cooling is very impor- 
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Fig. 6. — Test of the impact of HD cooling. The top panel shows 
the results using primordial abundances while the bottom shows 
for abundances from the cooling test. Initially, the temperature is 
started at 500 K and evolved for 100 Myrs. 

tant in both cases. 

3. MODEL FRAMEWORK 

Having developed and tested the chemistry and cool- 
ing routines necessary to study minihalo-shock interac- 
tions, we next turn to the detailed shock- minihalo inter- 
actions. Here we restrict our attention to a Cold Dark 
Matter (CDM) cosmology, with parameters are h = 0.7, 
Qo = 0.3 , Qa = 0.7, and ^5 = 0.045 {e.g. Spergel et 
al. 2007), where h is the Hubble constant in units of 100 
km s~^ Mpc~^, Qq^ r^A,and (^5 are the total matter, vac- 
uum, and baryonic densities, respectively, in units of the 
critical density. For our choice of /i, the critical density 

is Peru = 9.2 X 10-30 g/^^3^ 

3.1. The Minihalo 

A simple model is used for the gas and dark matter 
of the protocluster whose collapse redshift of Zc = 10 
(a cosmic age of ~ 0.5 Gyr) is taken to be just before 
the epoch of reionization, and whose total mass of Mc = 
3.0 X IO^Mq is taken to be on the large end of minihalos 
formed at this redshift. The gas is assumed to have a pri- 
mordial composition of 76% neutral atomic hydrogen and 
24% neutral atomic helium by mass. Initially, the cluster 
has a mean density that is enhanced by a factor A = 178 
(e.g. Eke, Navarro, & Frenk 1998) above the background, 
Pc = Al]o(l + Zc)^pcrit = 6.54 X 10"^^gm/cm3. In this 
case, the cloud's virial radius is Re = 0.393 kpc and virial 
velocity of ^c = 6.55 km/s. We assume that the radial 
profile is given by Navarro et al. (1997) 



p{R) 



VLqPc 



cx(l^cxY 3F{c) 



gm/cm^ 



(16) 



where c is the halo concentration factor, x = R/ Re , 
and Fit) = ln(l -\- 1) — j^. We assume that as the gas 
collapses inside the dark matter halo, it is shock- heated 
to its virial temperature, Tc = 1650 K and develops a 



density distribution of isothermal matter in the CDM 
potential well: 



Pgas(^) = Poe 



gm/cm^. (17) 



where the escape velocity is as a function of radius is 
given by v{xR^ir) = 2v'^[F{cx) -\-cx{l -\- cx)~'^][xF{c)]~'^ . 
From Madau et al. (2001), we take a typical value of 
the halo concentration to be c = 4.8, although some ob- 
servations suggest that high-redshift haloes maybe less 
concentrated than expected from this estimate (Bullock 
et al. 2001). With this value of c, we can compute the 
central density as: 



Po 



(178/3)0^1766^ (1 + zY 



j^ii^tru^dt 

39215 l^^Pcrit (1 + Zcf 

2.16 X 10"^^gm/cm^, 



-Pcrit 



(18) 



where A = 2c/F{c) = 10.3 and t = ex. 

To determine which case to use in our chemistry rou- 
tine the optical depth for H+, He+, and He++ recombi- 
nation was calculated from this profile: 



ru{r) = / Gj,n{r')dr' , 

Jrn 



(19) 



where a^, is the cross section of interaction and n{r^) is 
the number density. For hydrogen-like atoms the cross 
section is 



7.91 X 10-1^ 

Z2 



m 



g cm 



(20) 



where hui = 13.6 Z^ eV and Z is the nuclear charge. 
Calculating the optical depth from eqs. (17 - [20| yields 
the results given in Table 1. We assume fortne case 
of He++ recombination, that the surrounding helium is 
singly ionized. In all cases the optical depth is much 
greater than 1 (see §3.0) and therefore we use case B 
rates for all recombination reactions in our fiducial sim- 
ulations. 



TABLE 1 

Optical Depths, r is the optical depth to the center of 

THE CLOUD. In all CASES THE OPTICAL DEPTH IS MUCH GREATER 
THAN 1. 



Species 


Incident Energy eV 


r 


H+ 
He+ 
He++ 


13.6 
24.6 
54.4 


3553.9 
1049.2 
280.6 



3.2. Gravity 

As the minihalo in our simulation is made up of both 
gas and collisionless dark matter, a two-part gravity 
scheme was required. First, the Ricker (2008) multi- 
grid Poisson solver was used to calculate the gravita- 
tional potential due to the gas component. Second, 
the acceleration due to the total matter was calculated 



using Eq. (17) above. The general equation for the 
gravitationalacceleration of an ideal gas is given by 



^ [nR) 



dp{R) 
OR 



-p{R) 



dTjR) 
OR 



cm/s^ . To ac- 



<^Grav = pi^R)mp 

count for the dark matter halo's contribution to gravity, 
we calculated the gravitational acceleration of the total 
matter and subtracted its contribution from the bary- 
onic matter in the initial configuration using the above 
equation with constant temperature. Finally, we added 
the gravitational contribution from the self-gravity of the 
gas. The total acceleration is given simply by: 



<^Tot = <^M,0 — <^gas,0 + CLSG, 



(21) 



where aM is the acceleration from the total initial mass 
density as given by Eq. (17), agas,o is the contribution 
from the baryonic matter in the initial configuration, and 
asG is the contribution from the self-gravity calculated 
from the Poisson solver. Initially, when the minihalo 
gas is in hydrostatic balance with its surroundings, these 
last two terms will cancel each other and the cloud will 
remain unchanged. When the cloud is disrupted and 
cooling takes effect, the self- gravity will cause the cloud 
to collapse. 

The above equations are correct up to the viral radius 
of the cloud. To ensure a smooth density transition from 
the cloud, we simply keep the gas outside of the cloud 
be gravitationally bound to the cloud and solve for the 
expected density profile. The acceleration here is then. 



ClGv. 



and the density is 



p{R > Re) 



GM{R > Re) 



p{Rc)e(^ 



(22) 



(23) 



with Rq = GMcfrip/kbT. Here G is the gravitational 
constant, mp is the mass of a proton, ki, is Boltzmann's 
constant, and, as above Mc and Tc are the mass and tem- 
perature of the cloud. As R ^ oo, the density goes down 
to a small fraction of p{Rc)- A test of our gravity routine 
showed that the cloud was able to maintain hydrostatic 
balance for many dynamical times in the absence of an 
impinging galaxy outflow. 

3.3. The Outflow 

A Sedov- Taylor solution is used to estimate the prop- 
erties of the galactic outflow. The initial input energy 
is taken to be £^ = e£^55(ergs), where £^55 is the energy 
of the supernovae driving the wind in units of 10^^ ergs, 
and the wind efliciency e is derived from the amount of 
kinetic energy from the supernovae that is channeled into 
the outflow. The shock expands into a gas that is S times 
greater than the background at a redshift of Zc. As in 
Scannapieco et al. (2004), we assume that the cloud is 
a distance Rs = 3.6 kpc, using flducial values: Zc = 10, 
^55 = 10, Me = 3, 5 = 44 and e = 0.3 

With these values, the velocity of the blast front is 
Vs = 415 km s~^, when it reaches the minihalo, and 
the resulting temperature of the fully-ionized post-shock 
medium is T = 2.4 x 10^ K. By the time the shock 
has covered the separation distance, Rs, it will have en- 
trained a mass 



with a surface density of 

as = 2.6 X 10^ Mo kpc-^ 



(25) 



Note that while the above equations are the solutions 
that come from a simple spherical blast wave, the wind 
in the simulation is still well-approximated by a plane 
wave solution, because the size of cloud is much smaller 
than the distance between the supernova and the cloud. 
To model this wind, a time-dependent boundary con- 
dition is imposed at the leftmost boundary of the sim- 
ulation volume. The expected lifetime of the shock is 
given by ds = Vpost Ppost ts, where i;post is the post-shock 
velocity of the blast wave, Ppost is the post-shock density, 
and as is the surface density of the entrained material. 
Solving for ts and putting in the appropriate values, the 
expected shock life time is ts = 2.5 Myr. After this time 
the shock begins to taper off with the density decreas- 
ing and temperature increasing and keeping the pressure 
constant. This is done to prevent the excessive reflne- 
ment that a sharp cutoff would cause. The density falls 
off as 

Pit) _r.n. , nan.-r,/1.5 (26) 



P(0) 



= 0.01+0.99e-^'/^-^ 



and the temperature rises as 



no) 

T{t) 



= 0.01 + 0.99e-^^/i-^ 



(27) 



M,,Total = 4.4 X 10^ Mo, 



(24) 



where r^ is deflned as ^ ^y^ Y ' This also prevents the 
hydrodynamic (Courant) timestep from becoming ex- 
tremely short behind the shock, in order to maintain 
pressure equilibrium in an extremely rarifled medium. 

Note that in this initial study, the dark matter and gas 
distributions have been somewhat idealized, and more 
complicated geometries could be used to model these 
components in greater detail. For example, a triaxial 
instead of a spherical distribution could be assumed for 
the dark matter halo, inhomogeneities could be added 
to the minihalo gas, and the shock could be assumed to 
impact the minihalo off-axis. While each of these possi- 
bilities would be qualitatively interesting, and naturally 
alter the flnal outcome of the halo, they are nevertheless 
beyond the scope of this study. 

4. RESULTS 

Our simulations were carried out in a rectangular box 
with an effective volume of 3.2 xlO^ pc^. The ?/-axis and 
z-axis were the same length of 1170 pc and range be- 
tween [-585, 585] pc while the x-axis was twice as long, 
stretching between [-585, 1170] pc. The shock started on 
the left boundary while the cloud was centered at [0,0,0] 
pc. As hydrodynamic reflnement criteria, FLASH uses 
the second derivative of "reflnement variables," normal- 
ized by their average gradient over a cell. If this was 
greater than 0.8, the cell was marked for reflnement, and 
if all the cells in a region lie below 0.2, those cells were 
marked for dereflnement. 

A detailed summary of the runs performed is given in 
Table 2. The runs are labeled as either high or low reso- 
lution (H or L), whether atomic H-He recombination fol- 
lows Case A or Case B (A or B), and whether we impose 
a UV background (Y or N). The high-resolution. Case B, 
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no-background run (HBN) is taken to as our fiducial run 
and compared against other choices of parameters below. 



TABLE 2 

Summary of the numerical simulations in this study. 



Name 


Iref 


Resolution(pc) 


Cooling Mode 


Background (J21) 


HBN 


6 


4.55 


Case B 





LBN 


5 


9.11 


Case B 





HBY 


6 


4.55 


Case B 


10-1 


LBY 


5 


9.11 


Case B 


10-1 


HAN 


6 


4.55 


Case A 





LAN 


5 


9.11 


Case A 






4.1. Hydrodynamic Evolution 

In this simulation, several distinct stages of evolution 
are identified during the interaction between the cloud 
and the outfiow, as shown in Figures [7] and [SJ Initially, 
the cloud is in hydrostatic equilibrium as the shock en- 
ters the simulation domain. If it were not supported by 
pressure, the cloud would collapse on the free-fall time 
which, using the average cloud density, is 



te 



37r 
32Gp 



100 Myr. 



(28) 



As the cloud is initially in hydrostatic balance, the initial 
sound crossing time is similar to the free-fall time. 

As the shock contacts and surrounds the cloud, it heats 
and begins to ionize the gas. The shock completely en- 
velops the cloud on a characteristic "intercloud" crossing 
time scale, defined by Klein et al. (1994) as 



tu 



2Rc 



4.5 Myr. 



(29) 



As the cloud is enveloped, the shock moves through the 
outer regions fastest and ionizes this gas first. This in 
turn promotes rapid molecule formation as the gas cools 
and recombines incompletely, leaving 11+ and H~ to cat- 
alyze the formation of H2 and HD. Interestingly, because 
the shock slows down as it moves through denser mate- 
rial, the gas behind the center of the halo remains undis- 
turbed until the enveloping shocks meet along the axis at 
the back of the halo. The leads to a "hollow" H2 distri- 
bution at 6.6 Myrs, in which the molecular coolants are 
confined to a shell surrounding the undisturbed, purely 
atomic gas. 

After the enveloping shocks collide at the back of the 
cloud, a strong refiected shock is formed that moves away 
from the rear of the cloud and back through the halo 
material. Without cooling, this refiected shock would 
eventually lead to cloud disruption (Klein et al. 1994). 
However in our case, the shock has the opposite effect. 
It moves through the cloud, and the gas is briefiy ionized, 
but then quickly cools and recombines, forming H2 and 
HD throughout the cloud. This can be seen in the upper 
row of Figure [8J which shows the conditions at ~ 8 Myrs. 

At this point the cloud is denser, smaller, and full of 
new coolants. Using the conditions from the center of the 
cloud 8 Myr after the start of the simulation, we calculate 



new timescales. Now the freefall time is 21 Myr and the 
sound crossing time is ^ 27 Myr. The cloud is cold and 
dense enough to start collapsing. 

The timescale for the formation of H2, given in GA08, 
is 

*«■ = A^f* P°' 

where Xn^ is the mass fraction of H2, Xg is the mass 
fraction of electrons, ki is the reaction rate for the for- 
mation of H~ (H + e~ ^ H~ + 7), and n is the total 
number density (~ 1 cm~^ at 8 Myrs). Initially, as the 
shock begins to impact the cloud, this timescale is very 
short, on the order of 0.1 Myrs to get a final abundance 
of ~ 10~^. As the abundance of H2 increases and the 
abundance of electrons decrease, this timescale quickly 
increases. Although as the cloud collapses the density 
increases which lowers this timescale. 

The H2 cooling timescale, given by Klein et al. (1994) 
is 

l.bnkT 

^cool = 7 (S), (31) 

^H2^hAh,H2 

where nn^ and uh are the number densities of H2 and 
H respectively, and Ah,H2 ^^ ^^^ cooling rate between H 
and H2. At 8 Myr the H2 cooling time in most of the 
cloud is only 0.2 Myr, meaning that pressure support 
drops dramatically after this time. Any expansion due 
to shock heating is halted as the gas is quickly cooled 
by H2 and HD as they form. Furthermore, as the cloud 
collapses, the chemistry and cooling timescales decrease, 
rapidly accelerating the collapse. 

The final state of the cloud in our simulation is a thin 
cylinder stretching from the center of the dark matter 
halo to several times the initial virial radius. The tem- 
perature of this gas is 100 to 200 degrees, much colder 
than the initial virial temperature. The gas is also much 
denser than the initial minihalo, reaching values of up to 
10~^^ cm~^ or n ~ 10^ cm~^, in the center of the cloud, 
and even this density is probably only a lower limit set 
by the resolution of our simulation. On the other hand, 
the cloud is quite extended along the x-axis, with sub- 
stantial differences in velocity along the cylinder. Thus 
it is continually stretched and fragments until the end of 
the simulation at 14.7 Myrs (Row 3 in Figure [8|. 

Figure [9] shows rendered density contours of the major 
stages of evolution of the cloud from t = through the 
end of the simulation. The first panel shows the initial 
configuration, with the cloud in hydrostatic equilibrium, 
and the shock front entering from the left side of the 
simulation volume. The next panel shows the cloud af- 
ter being impacted by the shock, highlighting the density 
enhancement in the outer shell of the minihalo gas. The 
bottom left panel shows the cloud as it begins to cool 
and collapse, at a time at which the reverse shock has 
already passed though the cloud and coolants are found 
throughout the shocked, recombined material. Finally, 
the last panel shows the distribution at the end of the 
simulation. The cloud has now been stretched over a 
large distance and much of its mass has been acceler- 
ated to above the escape velocity, moving outside of the 
dark matter halo. The dense knots of this material in 
this figure are tightly gravitationally bound, have num- 
ber densities approaching 10^ cm~^, and are destined to 
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Fig. 7. — Initial evolution of the fiducial run, HBN, from t = through t = tic the time the shock completely surrounds the cloud. Each 
row shows the conditions in the central in a slice through the center of the simulation volume at times of (top), 2.1 (second row), 4.2 
(third row), and 6.6 Myrs (bottom row). The first column shows contours of the log of density from pgas = 10"-^^ to 10"^-*^ g cm"*^, which 
corresponds to number densities from n ^ 10"-^ to lO"^, the second column shows contours of the log of temperature from T = 10 to lO^K, 
and the third column shows contours of the log of the H2 mass fraction from Xn^ = 10~^ to 10"-*^. 



form extremely compact stellar clusters. 

4.2. Stellar Clusters 

While the collision happens on order of the shock cross- 
ing time of the halo, the final distribution of the clumps 
evolves on the longer timescale defined by i^vir/^c ~ 100 
Myrs. To study the final state of the stretched and col- 
lapsed distribution without continuing the simulation out 
to such extremely long times, we divided the x-axis into 
100 evenly spaced bins between x = kpc and x = lA 
kpc. We then calculated the mass of each bin by sum- 
ming up the gas from each cell from the FLASH simula- 
tion in a cylinder with a 24 pc radius and length of the 
bin. Similarly, we calculated the initial velocity of each 
bin by adding the momentum from each cell within this 
cylinder and dividing by the total mass in each bin. 

We evolved this distribution forward in time using a 
simple numerical model, which assumed that motions 
were purely along the x axis and pressure was negligi- 



ble at late times. In this case, acceleration could be cal- 
culated directly from the gravity between each pair of 
particles and from the potential of the dark matter halo. 
Furthermore, if any given particle moved past the parti- 
cle in front of it we merged them together, adding their 
masses and calculating a new velocity from momentum 
conservation. 

Evolving the distribution in this way for an additional 
200 Myrs past the end of the simulation yielded the re- 
sults shown in Figure [lOJ As the stretched cloud con- 
tinues to move outwaroTparticles begin to attract each 
other, and eventually merge together to create larger 
clumps. By 100 Myrs most of the particles have merged, 
after which their motions are purely ballistic. This can 
seen in the top panel as the lines for the late times over- 
lap each other and in the middle panel as the velocity 
profiles overlap. 

At the final time of 200 Myrs after the end of numeri- 
cal simulation, three small, stable clumps with masses of 
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Fig. 8. — Final evolution of the cloud from the propagation of the reverse shock across the cloud at t = 7.7 Myrs (top row), to collapse 
dit t = 11.8 Myrs (center row), through to the end of the simulation at t = 14.7 Myrs (bottom row). The panels have been cropped to show 
only the extended mass along the x-axis. Columns, values, and contours are the same as Fig. |7| 




Fig. 9. — Rendered density snapshots from the fiducial run, HBN. Colors show contours of log p from 10~^^ to 10"^-*^ cm~^. The top 
left panel, t = 0.0 Myrs, shows the initial setup with the stationary minihalo and the shock entering on the left. Top right, t = 6.3 Myrs, 
shows the state of the cloud as the shock as it envelopes the cloud. Bottom left, t = 9.4 Myrs, shows the cloud during collapse and cooling. 
Finally, bottom right, t = 14.7 Myrs, shows the final state of the cloud as it is stretched. Dense clumps can be seen in this panel, which 
we expect to become compact stellar clusters. 
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Fig. 10. — Evolution of the cloud up to 200 Myrs after the end of 
the simulation. The x-axis in each panel is the cumulative mass in 
solar masses. The top panel shows the mass of each particle, the 
middle panel shows their velocities, and the bottom panel shows 
their positions. The solid green lines show the profile at the end 
of the simulation tf = 14.7 Myrs, the dotted blue lines show the 
profile 50 Myrs later, the short-dashed cyan lines show the profile 
at tf + 100 Myrs, the dot-short dashed magenta lines show the 
profile at tf + 150 Myrs, and finally the short dash-long dashed red 
lines show the profile at tf + 200 Myrs. As time progresses we find 
that much of the material in the linear feature from Fig. [s] merges 
together. Most of this merging is complete by 100 Myrs after the 
end of the simulation. 

5.0x10"^ Mo, 4.0x10"^ Mo, and 3 x 10"^ Mo, as can be 
seen in the top panel of Fig. [lOJ Each of these new peaks 
is located far outside of the original dark matter halo. 

4.3. Case A vs. Case B 

At temperatures above 10^ K, the primary source of 
cooling is atomic lines from hydrogen and helium. Al- 
though we have shown that for the primordial cloud. 
Case B rates should be used for both the chemical net- 
work and cooling functions. Figure [TT] shows a compari- 
son between our fiducial run, HBN, and a run in which 
reaction and cooling rates are taken for case A recombi- 
nation (HAN). The high temperature Hydrogen- Helium 
cooling curve is taken from Weirsma et al. (2009). The 
upper panels show density contours and the bottom show 
contours of H2 abundance. 

As expected, the Case B simulation produces greater 
molecular coolant abundance at similar overall densities. 
This difference can be seen in the lower two panels of 
the first column of Fig. [TT] at 6.63 Myrs. To remain in 
pressure support as the cloud becomes denser from the 
shock, the cloud must get hotter. However, because Case 
A cools slightly faster, this support is quickly removed 
and the cloud takes on a more extended shape as evi- 
dent in the first two columns of Fig. [TT] Although, by 
14 Myrs the abundance of molecular coolants are very 
similar between HBN and HAN with each containing 



In both cases, the fate of the minihalo gas is the 
same. Atomic cooling occurs sufficiently rapidly to sap 
the shock of its energy and drop the post-shock temper- 
ature to ^ 10^ K, and nonequillibrium processes step 
in to provide molecular coolants below 10^ K. The gas 
is then able to collapse and form into a long dense fila- 
ment within which clumps are formed. In fact the only 
substantial differences between the runs are the details 
of the distributions of clumps, which is somewhat more 
extended in the case A run as compared to the case B 
run. 

4.4. UV Background 

A more uncertain aspect of our simulation is our as- 
sumption of a negligible dissociating background. In fact, 
the presence of at least a low level of dissociating back- 
ground is necessary in order for the minihalo not to col- 
lapse and form stars on its own, cooling by H2 and HD 
left over from recombination. To set an upper limit on 
the impact of such a background we modify the rates in 
our chemical network to approximate a relatively large 
dissociating background of J21 = 0.1, as discussed in sec- 
tion 2.1. Furthermore, as these rates are modified for all 
reactions throughout the simulation, this background is 
taken to affect even the densest regions of the cluster. 
This is equivalent to assuming that the cloud is optically 
thin to 11.2 to 13.6 eV photons at all times during the 
simulation. 

Figure [12] shows the comparison between the run in- 
cluding this background (HBY) and the fiducial run 
HBN. As expected, the abundance of H2 is reduced in 
the case with the UV background, peaking at about 
^ 10""^ instead of ~ 10~^ in the run without a back- 
ground. Interestingly, this difference persists even after 
a few megayears into the simulation, and the abundance 
of the HPY run remains stable at about ~ 10~^. 

However, this value is more than sufficient to cool the 
gas to the same temperature as in HBN. Even with this 
lower mass fraction, the cooling time is smaller than 
the dynamical time, and the evolution of the cloud re- 
mains essentially unchanged. The cloud collapses and is 
stretched into the same configuration as found without 
a background. Dense clouds are again found between 
0.2 and 0.4 kpc, at 0.55 kpc, and 0.9 kpc and the den- 
sity and temperature of each of these clouds is compara- 
ble to those found in the fiducial without a background. 
By neglecting any molecular self-shielding, this repre- 
sents a worst case scenario for H2, yet shock minihalo- 
inter actions continue to make compact stellar clusters. 

4.5. Resolution 

Finally, we considered the impact of the maximum re- 
finement level on our results. Figure 13 compares the 
fiducial run with 6 levels of refinement and an effective 
resolution of 4.55 pc to a lower resolution run (LBN) with 
5 maximum levels of refinement and an effective resolu- 
tion of 9.1 pc. Here density is shown in the upper two 
rows in each pair while the mass fraction of H2 is shown 
in the bottom two rows. 

Only a slight dependence on the formation for H2 is 
found with resolution. This is most apparent in the sec- 
ond column, corresponding to t = 7.7 Myrs, which shows 
that the shocks is slightly broadened in the lower reso- 
lution case. Since both chemical reactions and cooling 



Gray & Scannapieco 



13 



1.860543 Myrs 



'3.953041 Myrs 




Fig. 11. — Comparison between Case A and Case B cooling and chemistry rates. In this plot time varies across columns, moving from 
t = 6.6 Myrs (left column), to 7.7 Myrs (center column), to 14.0 Myr (right column). The upper two rows show the density in the central 

to 
10~2i g cm~^. The lower two rows show the H2 mass fraction in the fiducial run (third row) and the Case A run (bottom row). Here the 



slice from the fiducial, case B run (HBN, top row), and the case A run (HAN, second row), with log contours ranging from p = 10~^^ 
-21 g cjji~^. The lower two rows show the t 
; H2 mass fraction contours range from Xn^ 



10-8 to 10-1. 

go as n^, this smearing out has the effect of shghtly de- 
creasing H2 formation and coohng in the lower resolution 
run. However, enough H2 is produced in both cases for 
the cloud to collapse efficiently, and evolve in the same 
manner up until late times, when the difference in H2 
abundance is small. 

Furthermore, we also conducted similar resolution 
studies using Case A recombination, and also modify- 
ing chemistry and cooling to account for the presence of 
a dissociating UV background. Again comparisons be- 
tween the high-resolution runs (HAN and HBY) with 
the low resolution runs (LAN and LBY) uncovered only 
weak differences with resolution. Compact stellar clus- 
ters were formed in all cases. 

4.6. Source of Halo Globular Clusters? 

If indeed shock-minihalo interactions lead to massive 
clusters of stars, the longest-lived stars in these objects 
will remain observable down to low redshift, providing 
a direct observational constraint on our results. Fur- 
thermore, the clusters formed in our simulations are ex- 
tremely compact, and thus unlikely to be tidally dis- 
rupted as they eventually become gravitationally bound 
to the ever larger structures that form over cosmological 
time. As discussed in Scannapieco et al. (2004), the most 
natural candidate for these old and dense stellar clusters 



is the population of metal-poor globular clusters, associ- 
ated with the halos of galaxies {e.g. Zinn 1985; Ashman 
& Bird 1993; Larsen et al. 2001; Strader et al. 2005; 
Brodie et al. 2006). 

There are several detailed properties of the clusters 
in our simulations that support this association. Halo 
globular clusters, like the more metal-rich (disk) glob- 
ular clusters, are extremely compact, with typical half- 
light radii of ~ 3 pc {e.g. Jordan et al. 2005). Unlike 
higher- met allicity globular clusters, however, which may 
have mostly formed at intermediate redshift s (Elmegreen 
2010) the age of all metal-poor globular clusters is be- 
tween 10 — 13 Gyrs, placing their formation within or 
before reionization, during the epoch in which minihalos 
had not yet been evaporated by ionization fronts. 

While globular clusters exits at a range of masses, 
their mass distribution is well described by a Gaussian 
in log;Lo(^^c) with a 0.5 dispersion, and mean at 10^ so- 
lar masses, precisely spanning the masses of the clumps 
formed in our simulations. Although the lower mass limit 
of globular clusters may be set by destruction through 
mechanical evaporation {e.g. Spitzer & Thuan 1972) and 
shocking that occurs as the cluster pass through the 
host galaxies {e.g. Ostriker et al. 1972), the upper mass 
limit appears to be an intrinsic property of the initial 
populations {e.g. Fall & Rees 1985; Peng & Weisheit 
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6.6275920 M- 



6739405 Myrs 





Fig. 12. — Comparison between the fiducial run (HBN) and a run including a dissociating background (HBY) at three important stages 
of evolution. As in Fig. |ll[ from left to right the columns correspond to t = 6.6, 7.7, and 14.0 Myrs. From top to bottom the rows represent 
log density contours in the fiducial run (Row 1) and the dissociating background run (Row 2), contours of log H2 mass fraction in run HBN 
(Row 3) and HBY (Row 4), and contours of log temperature from run HBC (Row 5) and HBY (Row 6). The limits of each panel are the 
same as Fig|7| The addition of a background greatly reduces H2 but has almost no effect on the dynamics of the interaction. 



1991; Elmegreen 2010). Furthermore, this upper hmit 
of ~ 10^ Mq in stars corresponds roughly with the 10^ K 
virial temperature above which atomic coohng becomes 
effective in high-redshift viriahzed clouds of gas and dark 
matter (Fah & Rees 1985). No ^ lO^M© stehar clusters 
can result from shock minihalo interactions, because no 
minihalos exist with gas masses ^ IO^Mq. 

The second clue as to the origin of globular clusters 
comes from observations of stars being tidally stripped 
from these objects. If, like galaxies, globular clusters 
are contained within individual dark matter halos, the 
stripping of stars would be highly suppressed, due to the 
increased gravitational potential. Surprisingly, no such 
evidence of dark matter halos is seen (Moore 1995), sug- 
gesting that these clusters formed through a mechanism 



markedly different than star formation within galaxies. 

Again, the triggered-star formation seen in our shock- 
minihalo simulations provides a natural explanation of 
this key observed property. While gas is initially gath- 
ered together by a collapsed dark matter potential, the 
lack of effective coolants prevents this collapse from con- 
tinuing to the point that stars are formed at the center 
of the potential. Instead the momentum of the imping- 
ing galaxy wind is such that it is able to accelerate the 
gas above the escape velocity of the halo - even as it in- 
duces cooling and compresses the gas to the point that 
it remains gravitationally bound on its own. The trig- 
ger for star formation and the mechanism that removes 
the dark mater halo are one and the same, resulting in 
a population of stellar clusters that form at the moment 
they break free from their dark-matter enclosures. 
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Fig. 13. — Impact of maximum levels of refinement. Row 1 and Row 3 show the density and H2 contour s re spectively for the fiducial 
HBN while Row 2 and Row 4 show contours of density and H2 for LBN. Contour levels are the same as Fig. |11| Time is given at the top 
of each panel, and proceeds from t = 6.6 Myrs (left column) to t = 7.7 Myrs (center column) to 14.7 Myrs (right column). 



5. CONCLUSIONS 

Cosmological minihalos provide the initial building 
blocks out of which larger dark mater halos and galaxies 
form. During reionization, these minihalos surrounded 
all young galaxies and consisted of metal-free neutral 
atomic gas. The lack of molecules in these objects stalled 
their evolution at virial temperature below ^ 10^ K as 
cooling from atomic hydrogen becomes inefficient and 
any molecular coolants are dissociated by UV Lyman- 
Werner photons, too hot for star formation in such small 
objects. Only interactions that cause non-equilibrium 
molecule formation could allow the gas in minihalos to 
collapse and form stars. 

Previous work has used ionization fronts to create these 
conditions (Cen 2001), however, 3D hydrodynamic sim- 
ulations show that for either a stellar or quasar source, 
instead of creating clouds of H2, the cloud is instead com- 
pletely photo-evaporated (Iliev et al. 2005; Shapiro et 
al. 2006). However galactic outflows are another option 
for triggering star formation. Shocks provide the condi- 
tions for non-equilibrium chemistry without completely 
destroying the cloud. 

To simulate these interactions, we have added a pri- 
mordial chemistry network and associated cooling rou- 
tines into FLASH3.1. This network traces collisional 
ionization and recombination of hydrogen and helium as 



well as the formation of two primary coolants in the ab- 
sence of metals: H2 and HD. We have also included the 
impact of a dissociating background on these rates and 
implemented routines that control the cooling of gas in 
the presence atomic and molecular line cooling. 

With these code developments in place, we have been 
able to simulate the interaction between a galactic out- 
flow and a primordial minihalo that results in structures 
that are identified as proto-globular clusters. The shock 
fills two very important roles. First, it ionizes the neutral 
gas found in the minihalo, which recombines and begins 
to form H2 and HD, through nonequillbrium processes. 
These coolants allow the cloud to cool to much lower 
temperatures, triggering star formation. Secondly, the 
shock imparts momentum into the gas and accelerates it 
above the escape velocity. This creates a cloud of dense, 
cold molecular gas that is free from dark matter halos. 

Together, these many similarities suggest that there 
is not only room at low redshifts for the kinds of clus- 
ters were see in our simulations, but that low-redshift 
observations require the formation of stars by a mech- 
anism remarkably similar to the one we are simulat- 
ing. This physical connection requires further investi- 
gation, especially in light of a final remarkable property 
of globular clusters, the chemical homogeneity of their 
stars. While there are large metallicity between differ- 



16 



ent globular clusters the dispersion of [Fe/H] within any 
given cluster is less than 0.1 dex, or a factor of ^ 1.25 
(Suntzeff 1993). Usually, this is explained either by pre- 
enrichment, meaning that stars form out of gas that had 
already been homogeneously enriched by a previous gen- 
eration of supernovae {e.g. Elmegreen & Efremov 1997; 
Bromm & Clarke 2002), or "self-enrichment," meaning 
that the protocluster cloud was enriched by one or more 
SN events occurring within it (e.g.. Brown et al. 1995; 
Cen 2001; Nakasato et al. 2000; Beasley et al. 2003; Li 
& Burstein 2003). Both types of scenarios have strong 
flaws. In the pre-enrichment case, the key problem is 
just exactly what the previous population was, why it 
played only a secondary role in the formation history of 
the cluster, and how it could have enriched this material 
on very short timescales. In the self-enrichment picture, 
on the other hand, the key problems is that the kinetic 
energy corresponding to resulting SN is likely to unbind 
the gaseous protocluster (Peng & Weisheit 1991) long be- 
fore it enriches the gas. Recent two dimensional numeri- 
cal simulations of primordial supernova inside minihalos 
show that they do in fact unbind these systems over a 
range of supernova types and halo sizes (Whalen et al. 
2008b). 

Shock minihalo interactions have the potential to com- 
bine the best features of both these scenarios, bringing 



in the metals from an outside source, such that the gas 
collapses instead of unbinds, but nevertheless explaining 
how star formation could be synchronized in a large mass 
of coolant-filled gas. Yet highly homogenous population 
of stars can only be formed if the incoming, enriched gas 
and the primordial, minihalo material are sufficiently well 
mixed before star formation begins. Modeling this pro- 
cess will require simulations that build on the ones here 
to include metal line cooling above and below 10^ K and 
detailed models of subgrid-mixing. Testing this assump- 
tion and determining what sets of parameters lead to 
realistic globular cluster formation will be the focus of 
the upcoming papers in this series. 



We are grateful to Tom Abel, Jean Brodie, Simon 
Glover, Raul Jimenez, Joaquin Prieto, Jay Strader, Sum- 
ner Starrfield, F. X. Timmes, and Steven Zepf for their 
helpful comments. ES acknowledges the support from 
NASA theory grant NNX09AD106. Ah simulations were 
conducted on the "Saguaro" cluster operated by the 
Fulton School of Engineering at Arizona State Univer- 
sity. The results presented here were produced using 
the FLASH code, a product of the DOE ASC/Alhances- 
funded Center for Astrophysical Thermonuclear Flashes 
at the University of Chicago. 



REFERENCES 



Abel, T., et al 2002, Science, 295, 93 

Ahn, K., Shapiro, P.R., Iliev, I.T., Mellema, G., & Pen, U. 2009, 

ApJ, 695, 1430 
Ashman, K. M., &; Bird, C. M. 1993, AJ, 106, 2281 
Bader, G., &; Deuflhard, P 1983, Numer. Math. 41, 373 
Barkana R., &; Loeb A. 1999, ApJ, 523, 54 
Beasley, M. A., Kawata, D., Pearce, F. R., Forbes, D. A., & 

Gibson, B. K. 2003, ApJ, 596, L187 
Bond J. R., Szalay A. S., &; Silk J. 1988, ApJ, 324, 627 
Brodie, J.P., &; Strader, J. 2006, ARA&A, 44, 193 
Bromm, V., Coppi, P.S., &; Larson, R. 2002, ApJ, 564, 23 
Bromm, V., & Clarke, C. J. 2002, ApJ, 566, LI 
Brown, J. H., Burkert, A., &; Truran, J. W. 1995, ApJ, 440, 666 
Bullock J.S, et al. 2001, MNRAS, 321, 559 
Cen, R. 2001, ApJ, 560, 592 

Ciardi, B., & Ferrara, A. 2005, Space Sci. Rev., 116, 625 
Ciardi, B., Ferrara, A., &; Abel, T. 2000, ApJ, 533, 594 
Colella, p., & Glaz, H.M. 1985, JCoPh, 59, 264 
Colella, p., &; Woodward P. 1984, JCoPh, 54, 174 
Dalgarno, A., &; McCray, R. A. 1972, ARA&A, 10, 375 
Eke, v., Navarro, J., & Frenk, C. 1998, ApJ, 503, 569 
Elmegreen, B. G., 2010, in lAU Symp. 226, Star Clusters Basic 

Galactic Building Blocks Throughout Time and Space, ed. R. 

de Grijs &; J. Lepine, Cambridge, 3 
Elmegreen, B. G., & Efremov, Y. N. 1997, ApJ, 480, 235 
Fan, S.M., & Rees, M.J. 1985, ApJ, 298, 18 
Ferland, G.J., Korista, K.T., Verner, D.A., Ferguson, J.W., 

Kingdon, J.B., &; Verner, E.M. 1998, PASP, 110, 761 
Ferrara, A. 1998, ApJ, 499, L17 

Fragile, C.P., Murray, S.D., &; Anninos, P. 2004, ApJ, 604, 74 
Franx, M., Illingworth, G. D., Kelson, D. D., van Dokkum, P. G., 

& Tran, K.-V. 1997, ApJ, 486, L75 
Fryxeh, B., et al. 2000, ApJS, 131, 273 
Fryxeh, B., Miiller, E., &; Arnett. B. 1989, nuas.conf, 100 
Fujita, A., Martin, C. L., Mac Low, M.-M., & Abel, T. 2003, 

ApJ, 599, 50 
Gain, D., & Palla, F. 1998, A&A, 335, 403 
Glover, S.C.O., & Abel, T. 2008, MNRAS, 388, 4 (GA08) 
Glover, S.C.O. 2009, private communication (G09) 
Glover, S.C.O. , Savin, D.W. 2009, MNRAS, 393, 991 
Haiman, Z., Able, T., &; Madau, P. 2001, ApJ, 551, 599 
Haiman, Z., Rees, M., &; Loeb, A. 1997, ApJ, 476, 458 (erratum 

484, 985) 



Heckman, T. M., Lehnert, M. D., Strickland, D. K., & Armus, L. 

2000, ApJ, 129, 493 
Iliev, L T., Shapiro, P. R., &; Raga, A. C. 2005, MNRAS, 361, 405 
Jordan, A. et al. 2005, ApJ, 634, 1002 
Kang, H., Shapiro, P.R., Fah, M., &; Rees, M.J. 1990, ApJ, 363, 

488 
Kaps, P., & Rentrop, P. 1979, Numer. Math., 33, 55 
Klamer, L J., Ekers, R. D., Sadler, E. M., & Hunstead, R. W. 

2004, ApJ, 612:L97-L100 
Klein, R., McKee, C.F., &; Colella, P. 1994, ApJ, 420, 213 
Larsen, S. S., Brodie, J. P., Huchra, J. P., Forbes, D. A., &; 

Grillmair, C. J. 2001, AJ, 121, 2974 
Lehnert, M. D., &; Heckman, T. M. 1996, ApJ, 462, 651 
Li, Y., &; Burstein, D. 2003, ApJ, 598, L103 
Lipovka, A., Nunez-Lopez, R., &; Avila-Reese, V. 2005, MNRAS, 

361, 850 
Mac Low, M., &; Shuh, J.M. 1986, ApJ, 302, 585 
Madau, P., Ferrara, A., &; Rees, M. 2001, ApJ, 555, 92 
Martin, C. L. 1998, ApJ, 506, 222 
Martin, C.L. 1999, ApJ, 513, 156 
Moore, B. 1995, ApJ, 461, L13 

Nakasato, N., Mori, M., & Nomoto, K. 2000, ApJ, 535, 776 
Navarro, J.F., Frenk, C.S, & White, S.D.M. 1997, ApJ, 490, 493 
O'Shea, B.W., & Norman, M.L. 2007, ApJ, 654, 66 
Osterbrock, D.G., Astrophysics of Gaseous Nebulae and Active 

Galactic Nuclei. University Science Books. 1989. 
Ostriker, J. P., Spitzer, L., &; Chevalier, R. A. 1972, ApJ, 176, L51 
Palla, F., Salpeter, E.E., &; Stahler, S.W. 1983, ApJ, 271, 632 
Palla, F., &; Zinnecker, H. 1988, in lAU Symp. 126, The Harlow 

Shapley Symposium on Globular Cluster Systems in Galaxies, 

ed. J. E. Grindlay & A. G. D. Philip (Doredrecht: Reidel), 697 
Peng, W., &; Weisheit, J. C. 1991, PASP, 103, 891 
Pettini, M., Kellogg, M., Steidel, C. C, Dickinson, M., 

Adelberger, K. L., & Giavalisco, M. 1998, ApJ, 508, 539 
Prieto, J. P., Infante, L., &; Jimenez, R. 2009, arXiv:0809.276v2 
Ricker, P. M. 2008, ApJS, 176, 293 

Rupke, D. S., Veilleux, S, & Sanders, D. B. 2005, ApJS, 160,115 
Sokasian, A., Yoshida, N., Abel, T., Hernquist, L., &; Springel, V. 

2004, MNRAS. 350, 47 
Scannapieco E., Ferrara A., &; Madau P. 2002, ApJ, 574, 590 
Scannapieco, E. et al. 2004. ApJ, 615,29 
Shapiro, P.R., &; Kang, H. 1987, ApJ, 318, 32 



Gray & Scannapieco 



17 



Shapiro, P. R., Raga A. C, & Mellema G. 1997, in "Structure 

and Evolution of the Inter- galactic Medium from QSO 

Absorption Line Systems," p. 41 
Shapiro, P. R., Raga A. C, & Mellema G. 1998, Mem. Soc. 

Astron. Ital., 69, 463 
Shapiro, P. R., Iliev, I. T., &; Raga, A. C. 2004, MNRAS, 348, 753 
Shapiro, P. R., Iliev, I. T., Alvarex, M. A., Sz Scannapieco, E. 

2006, ApJ, 648, 992 
Smith, B., Sigurdsson, S., &; Abel, A. 2008, MNRAS, 385, 1443 
Spergel, D. N. et al. 2007, ApJ, 170, 377 
Spitzer, L., &; Thuan, T. X. 1972, ApJ, 175, 31 
Strader, J., Brodie, J. P., Cenarro, A. J., Beasley, M. A., & 

Forbes, D. A. 2005, AJ, 130, 1315 
Suntzeff, N., Mateo, M., et al. 1993, ApJ, 481, 208 



Timmes, F.X. 1999, ApJ, 124, 241 

Thacker R. J., Scannapieco E., Sz Davis M. 2002, ApJ, 202, 581 
Uehara, H. & Inutsuka, Shu-ichiro 2000, ApJ, 531, L91 
Veilleux, S., Cecil, G., &; Bland-Hawthorn, J. 2005, ARA&A, 43, 

769 
van Breugel, W., Filippenko, A.V., Heckman, T.,Miley, G. 1985, 

ApJ, 293, 83 
Weirsma, R.., Schaye, J., &; Smith, B.D. 2009, MNRAS, 393, 9 
Whalen, D., O'Shea, B.W., Smidt, J., &; Norman, M.L. 2008a, 

ApJ, 679, 925 
Whalen, D., van Veelen, B., O'Shea, B.W., & Norman, M.L. 

2008b, ApJ, 682, 49 
Wiita, Paul J. 2004, Ap&SS, 293, 235 
Zinn, R. 1985, ApJ, 293, 424 



